load surface_time;
S=length(x);
timehistint=linspace(timehist(1),timehist(end),length(timehist)); dt=timehistint(2)-timehistint(1);

lowcut=0;
surfacehisthf=surfacehist;
for s=1:S
    s
    %surfacehist=interp1(timehist,surfacehist,timehistint); surfacehist(1)=surfacehist(2);
    tempsurface=interp1(timehist,surfacehist(s,:),timehistint); tempsurface(1)=tempsurface(2);
    temptemp=interp1(timehist,avtemphist(s,:),timehistint);
    avtemphist(s,:)=bandpass(temptemp,lowcut,inf,dt);
    surfacehist(s,:)=bandpass(tempsurface,lowcut,inf,dt);
    %surfacehist(s,:)=surfacehist(s,:)-surfacehist(s,1);
    surfacehisthf(s,:)=bandpass(tempsurface,1/100,inf,dt);
    
    for d=1:length(depthrange)
        temptemp=interp1(timehist,temphist(s,:,d),timehistint);
        %temptemp=temptemp-mean(temptemp);
        temptemp=bandpass(temptemp,lowcut,inf,dt);
        temphist(s,:,d)=temptemp;
        if depthrange(d)<40e3
            temphist(s,:,d)=temptemp+nan;
        end
        
        tempP=interp1(timehist,Phist(s,:,d),timehistint);
        tempP=bandpass(tempP,lowcut,inf,dt);
        Phist(s,:,d)=tempP;
    end
end


[TIMEHIST X]=meshgrid(timehistint,x);

[Zgrid,Xgrid]=meshgrid(depthrange,x);

for time=1:length(timehistint)
    time
    
    currentsurface=surfacehist(:,time);
    figure(1); clf;
%     subplot(2,2,[1 3]); pcolor(X,TIMEHIST,-surfacehist); shading flat; xlabel('Distance/m'); ylabel('Time/Ma'); colorbar; hold on; plot([x(1) x(S)], [1 1]*timehistint(time),'k','linewidth',2);
%     subplot(2,2,2); plot(x,-surfacehist(:,time)); axis([0 2000e3 -170 170]); colorbar;
%     tempslice=reshape(temphist(:,time,:),S,length(depthrange));
%     subplot(2,2,4); pcolor(Xgrid,Zgrid,tempslice); shading flat; set(gca,'ydir','rev'); axis equal; caxis([-130 130]); colorbar; axis([0 2000000 0 710000]); %title([num2str(time/(60*60*24*365*1e6)) ' Ma ']);     

    subplot(3,2,[1 3 5]); pcolor(X,TIMEHIST,-surfacehist); shading flat; xlabel('Distance/m'); ylabel('Time/Ma'); colorbar; hold on; plot([x(1) x(S)], [1 1]*timehistint(time),'k','linewidth',2);
    subplot(3,2,2); plot(x,-surfacehist(:,time)); axis([0 2000e3 -170 170]); colorbar;
    tempslice=reshape(temphist(:,time,:),S,length(depthrange));
    Pslice=reshape(Phist(:,time,:),S,length(depthrange));
    subplot(3,2,4); pcolor(Xgrid,Zgrid,tempslice); shading flat; set(gca,'ydir','rev'); axis equal; caxis([-130 130]); colorbar; axis([0 2000000 0 710000]); %title([num2str(time/(60*60*24*365*1e6)) ' Ma ']);
    subplot(3,2,6); pcolor(Xgrid,Zgrid,Pslice); shading flat; set(gca,'ydir','rev'); axis equal; caxis(5e6*[-1 1]);  colorbar; axis([0 2000000 0 710000]); %title([num2str(time/(60*60*24*365*1e6)) ' Ma ']);
    
    print(gcf, '-djpeg100','-r100',['..\outfiles\surf_surf_temp' num2str(time)]);
end